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Kitaev’s compass model on the honeycomb lattice realizes a spin liquid whose emergent excita¬ 
tions are dispersive Majorana fermions and static Z 2 gauge fluxes. We discuss the proper selection 
of physical states for hnite-size simulations in the Majorana representation, based on a recent pa¬ 
per by Pedrocchi, Chesi, and Loss [Phys. Rev. B 84, 165414 (2011)]. Certain physical observables 
acquire large finite-size effects, in particular if the ground state is not fermion-free, which we prove 
to generally apply to the system in the gapless phase and with periodic boundary conditions. To 
illustrate our findings, we compute the static and dynamic spin susceptibilities for finite-size sys¬ 
tems. Specifically, we consider random-bond disorder (which preserves the solubility of the model), 
calculate the distribution of local flux gaps, and extract the NMR lineshape. We also predict a 
transition to a random-flux state with increasing disorder. 


I. INTRODUCTION 

Frustrated magnetism is an exciting field of research 
in condensed matter physics. Particular attention has 
been devoted to so-called spin-liquid states/^ In a strin¬ 
gent definition, these are zero-temperature states of local- 
moment systems with half-odd-integer spin per crystallo¬ 
graphic unit cell which are characterized by the absence 
of any spontaneous symmetry breaking. Typically, the 
low-energy description of such states involves non-trivial 
elementary excitations with fractional quantum numbers 
which are coupled to an emergent gauge field. 

In a seminal paper ,I^Kitaev proposed a model of quan¬ 
tum spins 1/2 on a two-dimensional honeycomb lattice, 
subject to a particular type of anisotropic exchange in¬ 
teractions, often dubbed “compass” interactions This 
model is exactly solvable, thanks to an infinite set of 
conserved quantities. It realizes a non-trivial spin-liquid 
state which, depending on the interaction parameters, 
can be gapless or gapped. Its elementary excitations are 
dispersing spinless “matter” fermions which are coupled 
to a frozen Z2 gauge field. , By now, many properties 
of the Kitaev model have been studied, including staticP 
and dynamic!^ spin correlations as well as the physics 
of isolated defects.l^lf^ In addition, variants of the Kitaev 
model on other lattices, both in twcPEH and thred^^lli^ 
space dimensions, have been discussed. In all cases, the 
most popular analytical treatment of the compass inter¬ 
actions utilizes a Majorana representation of spins. Sub¬ 
tleties in dealing with the cor respo nding enlarged Hilbert 
space have been pointed out.^ISliZl 

On the materials side, oxides of the family A2lr03, 
with magnetic iridium ions subject to strong spin-orbit 
coupling, have been proposed^^ to realize an exchange 
Hamiltonian of Kitaev type, supplemented by addi¬ 
tional spin-symmetric Heisenberg interactions. The re¬ 
sulting Heisenberg-Kitaev model has been investigated 
extensivelyl^sHUl While the spin liquid is stable to small 
admixtures of Heisenberg interactions, larger perturba¬ 
tions destroy it in favor of a variety of magnetically 


ordered phases. Experimentally, both Na2lr03 and 
Li2lr03 have been found to display magnetic order at 
low temperatures,l2SH22l and it has been speculated that 
pressure might be used to tune them towards the spin- 
liquid regime. However, the precise microscopic Hamil¬ 
tonian describing the magnetism in A2lr03 is under 
debate.l22l2M33] 

In this paper, we consider the honeycomb-lattice Ki¬ 
taev model with random-magnitude exchange interac¬ 
tions, i.e., bond randomness. The model remains exactly 
solvable and thus belongs to the rare cases of exactly solv¬ 
able random spin models in dimensions d> 2. (Brief dis¬ 
cussions of disorder in the Kitaev model have been given 
in Refs. [3 and [311 and a Kitaev-style chiral spin-liquid 
model with random exchange was considered in Ref. 1351 ) 
We shall utilize the Majorana-fermion representation to 
investigate the magnetic response of the bond-disordered 
Kitaev spin liquid, in particular the NMR lineshape. Dis¬ 
order is treated exactly via finite-size exact diagonaliza- 
tion. 

Particular attention is paid to the proper selection of 
physical states in the Majorana representation,li^ which 
results in a condition on the parity of matter fermion ex¬ 
citations. While this condition generally depends on both 
the flux configuration and the system geometry, we are 
able to prove that, for clean systems with periodic bound¬ 
ary conditions and interaction parameters in the gapless 
phase, this parity must always be odd in the flux-free 
sector. Hence, the physical ground state is not fermion- 
free, but contains one matter fermion excitation. As we 
will show, this implies large finite-size effects for many 
observables. As an aside, we point out that the ground 
state of the clean Kitaev model for certain small systems 
is not in the flux-free sector of the Z2 gauge field. For 
large systems, we predict a quantum phase transition, 
upon increasing bond randomness, from a flux-free to a 
random-flux ground state. 

The body of the paper is organized as follows: In Sec¬ 
tion^ we introduce the random-bond Kitaev model to¬ 
gether with its Majorana representation and the numer¬ 
ical solution in terms of free canonical fermions. The 
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required projection to the physical Hilbert space is sub¬ 
ject of Section III Section IV outlines the numerical cal¬ 
culation of the susceptibility. In Section |V] we briefly 
show numerical results for observables in the clean sys¬ 
tem, with focus on their finite-size behavior. General 
aspects of quenched bond disorder in the Kitaev model 
are discussed in Section lYll while concrete numerical re¬ 
sults are presented in Section |VII| The transition to the 
random-flux state is discussed in Section IVIIII A sum¬ 
mary closes the paper. Technical aspects of the physical- 
state selection are relegated to the appendix, as is the 
comparison of the Majorana and exact solutions for a 
small system of four unit cells. 


II. MODEL AND MAJORANA 
REPRESENTATION 

A. Random-bond Kitaev model 

The Kitaev modeP describes spin-1/2 degrees of free¬ 
dom at sites j of a honeycomb lattice which interact 
via Ising-like nearest-neighbor exchange interactions J“. 
The anisotropy direction in spin space, a = z, is 
coupled to the bond direction in real space, reflecting 
a strong spin anisotropy from spin-orbit coupling. We 
generalize the model to spatially varying, i.e., random, 
couplings, such that the Hamiltonian reads 

= - E - E - E (1) 

(b>a (b)^ 

where (t“ are Pauli matrices, and denotes an a = 
x,y,z bond as in Fig. In the clean case = J^, 
Jij = ^ Jfj = For isotropic couplings, = jv = 

= J, the model possesses a symmetry of combined 

real-space and spin rotations. 

In our simulations of bond disorder, the exchange cou¬ 
plings J“- will be drawn from uncorrelated box distribu¬ 
tions with mean value > 0, £ [J“ — A“, J“ + A“]. 

In a possible experimental realization in an insulating 
solid, disorder in the Jij arises from random lattice dis¬ 
tortions and/or chemical disorder on non-magnetic sites, 
both of which locally modify individual exchange paths. 


B. Majorana representation 

Following Kitaev’s solution,!^ we introduce four (real) 
Majorana fermions b^, , b^ and c. Defining df = ibfci, 

the original Hamiltonian in Eq. Q can be mapped to 

l~iu — ^ ^ ^ ^ (2) 

{iJ) 

where Uij = ib^'^b^'\ Uij = —Uji, and the summation is 
over all nearest-neighbor bonds. We follow the conven¬ 
tion that, when specifying Uij, i is located on sublattice 



FIG. 1: Honeycomb lattice with basis vectors ei ,2 and an 
illustration of the periodic boundary conditions, characterized 
by the cluster size Li ,2 and the twist parameter M. The figure 
corresponds to Li — L 2 = 3 and M = 2. 


A. The operators Uij, with eigenvalues Uij = ±1, com¬ 
mute with each other and with the Hamiltonian "Hu, i-C., 
the {uij} are constants of motion. A given set {uij} re¬ 
duces the Hamiltonian to a bilinear in the c Majorana 
operators: 

= 2 (-M^ ili) ■ 

Here M is an A x A matrix with elements 
and ca(b) is a vector of A Majorana operators on the 
A{B) sublattice. Hence the problem takes the form of 
non-interacting Majorana fermions coupled to a static 
Z 2 gauge held. 

The eigenmodes of "Hu can be found via singular-value 
decomposition of M, M = USV'^, where U and V are 
A X A orthogonal matrices, and S' is an A x A diagonal 
matrix containing the non-negative singular values of M. 
We define new Majorana operators according to 

(5/,..., 6)v) = (cA.i, • ■ ■ ,ca,n)U , 

(4) 

(&'/, . . . , 6)(r) = (cs.l, ■ ■ ■ , Cb,n)V . 

We may combine the transformation matrices U and V 
into a matrix , 

^ (v 0) ’ 

which is equivalent to (5“ defined in Eq. (4) of Ref. [T7] 
after re-ordering of both rows and columns. 

For a given set of {uij} the Hamiltonian now has the 
form "Hu = * where > 0 are the singu¬ 

lar values of M. It is convenient to combine the Majorana 
operators 5', b” into canonical fermions according to 

am = ^(&L + *C)- (6) 
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This eventually gives 

N 

"Hu = ^ €m{‘2al^am - 1) (7) 

m—l 

with the ground-state energy Eq = Eigen¬ 

states of the Hamiltonian Q can thus be understood as 
a direct product of “gauge” (u) and “matter” (a) degrees 
of freedom. 


C. Boundary conditions 

To avoid edge effects, the analytical discussion as well 
as the numerical calculations will be performed for finite- 
size systems with periodic boundary conditions. We will 
comment on open boundary conditions in Section |III C| 
below. 

As in Ref. [IZl we will restrict our attention to “rectan¬ 
gular” clusters of size N = Li x L 2 unit cells, with 2N 
spins, but allow for a geometric “twist” characterized by 
an integer M when imposing periodicity. Here, the torus 
is defined through the basis vectors Liei and L 2 e 2 +Mei, 
see Fig. In the isotropic case this represents the most 
general set of periodic boundary conditions for rectangu¬ 
lar clusters. 


In the Majorana representation, the loop (or flux) op¬ 
erators W can be expressed through the bond variables 
Uij; the same holds for their eigenvalues. For instance, 
the plaquette fluxes take the form 

Wp = U2lU23U43U45Ue5UQi . ( 10 ) 

As a consequence of gauge invariance, the fermion spec¬ 
trum {cm} and the ground-state energy Eq depend on 
the Uij only through the values of the fluxes, {Wp} and 
1Fi,2. 

For a translation-invariant system of sufficiently large 
size, the ground state is located^ in the flux-free sector, 
corresponding to all IT = 1 (i.e. all Uij = 1). In this sec¬ 
tor, the excitation spectrum of the hopping Hamiltonian 
Hu can be found using a Fourier transformation. De¬ 
pending on the anisotropy of the couplings, the system is 
either gapped or gapless, with the latter case including 
the isotropic point, . Here, the low-energy 

part of the spectrum consists of two Dirac cones similar 
to graphene. It is worth noting that the ground state of 
certain small systems is not in the flux-free sector; this 
will be further discussed in Section IV HI 


III. PHYSICAL MANY-BODY STATES IN THE 
MAJORANA DESCRIPTION 


D. Flux degrees of freedom 

For every closed loop C of the lattice , the Kitaev model 
Q features a conserved quantity For a loop C 

containing L sites labeled {1,2, ...,T}, the corresponding 
operator is 

rir *“1,2 ;i“l,2 *“2,3 *“2,3 *“11,1 *“J1,1 /'0^ 

Wc — <^1 0^2 ^2 ^3 ■ ■ ■ ^1 (”) 

where aij = cc, y, or z corresponds to the type of the 
bond connecting sites i and j. The eigenvalues of the Wc 
are Wc = ±1, each corresponding to a Z 2 flux. It is con¬ 
venient to introduce loop operators for the flux through 
each elementary plaquette of the lattice, 

ITp = (9) 


The Majorana representation of spins 1/2 is overcom¬ 
plete: The total Hilbert space of Hq has 4^^ states, as 
compared to 2^^ states forming the Hilbert space of Hk- 
First, the 2^+^ physical flux sectors are represented by 
2^^ link variables Uij, such that different configurations 
of {uij} correspond to the same flux sector. Second, 
within each flux sector there are 2^ states of the c Majo¬ 
rana fermions, to be compared to 2-^“^ physical states. 
This implies that the possible fermion+flux states can 
be grouped into “physical” and “unphysical” states,l2E^ 
with not all physical fermion+flux states corresponding 
to different spin states of the model Q . 

A. Projection 


with 1,...,6 labelling the sites of the plaquette under 
consideration. For periodic boundary conditions, there 
are two additional (“topological”) loop operators Wi ^2 
that wrap around the torus in the direction of the unit 
vectors ei 2 and are related to the flux through the torus 
holes. 

A system with N unit cells and periodic boundary con¬ 
ditions is characterized by (N — 1) independent plaquette 
fluxes Wp, due to the constrainiP ITp = 1. Together 
with the torus fluxes lTi .2 the total number of flux de¬ 
grees of freedom is {N + 1). Given that the dimension of 
the physical Hilbert space of Hk is 2^^, this implies that 
each individual flux sector consists of 2^~^ many-body 
states. 


We first summarize the Majorana state projection as 
outlined in Refs. 121171 and then present our extended re¬ 
sults in the following subsections. 

An eigenstate of Hk, |C)) satisfies the condition 
Dj 1^) = 1^), where Dj = —iajaj&j = 1. Written in 
terms of Majorana operators, we have Dj = bjSjbjCj, 
with eigenvalues ±1. Now, a physical eigenstate must 
satisfy Dj |fy = + |fy for all j. One can therefore define 
a projection V to the physical subspace of the Majorana 
Hilbert space according to 


2N 

p=n 





( 11 ) 









4 


In this subspace the original spin Hamiltonian, "Hkj and 
Kitaev’s Majorana Hamiltonian for the honeycomb lat¬ 
tice, "Hu, are equivalent. The operator Dj can be thought 
of as an Ising gauge transformation. Since the spin op¬ 
erators are gauge-invariant, their matrix elements in any 
gauge-fixed sector are identical to that in the physical 
gauge-invariant subspace .1^ 


The effect of the projection (111 to annihilate an un¬ 
physical state is easily seen by rewriting it aJ^^ 


V = S 




= SV^ 


0) 


( 12 ) 


where S symmetrizes over all gauge-equivalent subspaces 
while Vq projects out unphysical states. 

The operator D — J([^- Dj can be expressed in the Ma¬ 
jorana representation. After re-ordering the fermion op¬ 
erators - see Appendix B of Ref. |T7] - it can be brought 
into the form: 

b= (-i)"n^^ n as) 

3 {ij)^ (ij) 

Here, tTc = parity of the c (matter) Majo¬ 

rana fermions, and we followed the convention that sites 
labeled with odd (even) numbers belong to the A{B) sub¬ 
lattice (this differs from Ref. [T7)l . The exponent 6* is a 
consequence of the anticommutation relation of the Ma¬ 
jorana fermions and depends on the lattice geometry. For 
the boundary conditions in Fig. it read^^ 

e = Li+ L 2 + M{Li- M). (14) 


An alternative representation of the Majorana states 
uses local complex fermions. For each unit cell r one can 
construct one complex matter fermion 

/r = ^ [cA,r - iCB,r] (15) 


and three complex gauge fermions defined on the bonds 
emanating from site i on sublattice A: 


X 


cx. 

r 




(16) 


Then we have icA,rCB,r = 1 ~ 2/1/r such that we can ex¬ 
press the parity Tr^ as = (—1)'^^' with Nf = /l/r- 

Similarly, ibfb^j = iiij = 1 — 2(x“)lx? which yields 
n(ij) = (“1)^’^- This allows one to rewrite the oper¬ 


ator D (13l using the fermion numbers Nf and N^: 


b = (17) 


The condition for a state being physical, b = 2Vo — 1 = 
1, selects states with either even or odd total fermion 
number, depending on the geometry factor (—1)®. For 
fixed {uij} this eliminates half of the many-body states 
from the Hilbert space of 'Hu, as anticipated, and implies 


that fermions can only be excited pairwise. We note that 
the factor (—1)®, derived in Ref. [T71 does not seem to 
appear in earlier works.t^ 

To convert Eq. 0 into a more useful form, it is im¬ 
portant to distinguish the parity rfc of the c fermions 
from the parity tt = n::(i — 2al^am) of the eigenmodes 
Q>m §. Given that the c and a fermions are related via 
the canonical transformation Q", Eq. © , one finddl^ 

-ifc = det(Q“)7r . (18) 

Combining Eqs. ( [T7| and 0 the operator b reads 

£) = (-l)®det(Q“)(-l)^“(-l)^^ (19) 

with Na = X)m being the number of matter 

fermion excitations. 


B. Fermion parity for periodic boundary conditions 


In general, the value of D (191 depends in a combined 


fashion on the flux configuration, the boundary condi¬ 
tions, and the distribution of the coupling constants; con¬ 
crete examples were given in Ref. 1171 

Here we go one step further: For a translation-invariant 
system in the gapless phase, we are able to prove that in 
the flux-free sector we have (—1)^ det(Q“) = —1 indepen¬ 
dent on the system geometry. Details of this proof are 
given in Appendix Since the flux-free sector is char¬ 
acterized by = 0, the condition Z) = 1 translates into 

TT = (—1)^“ = —1, i.e., all physical states in the flux-free 
sector must have an odd number of a fermion excitations. 
Hence, the naive fermion-free state is not a physical state. 
This has consequences for the calculation of observables, 
as will be discussed below. 

On general grounds, we expect that a single fermion in 
an extended system of size N can cause only 1 /N effects 
on observables. Hence, the proper selection of physical 
states discussed here, albeit important for finite-size sys¬ 
tems, is not expected to influence typical observables in 
the thermodynamic limit. Indeed, in our calculations 
we find strong differences in the finite-size behavior of 
observables calculated with either physical or unphysi¬ 
cal states, but these differences diminish with increasing 
system size. However, for observables where 1/N cor¬ 
rections are crucial - this applies to quantum impurity 
problems - the state of affairs might be different; this will 
be investigated in future work. 


C. Fermion parity for open boundary conditions: 
dangling gauge fermions 

The considerations in Ref. im and the present section 
show that, for a Kitaev model with periodic boundary 
conditions, half of the Majorana many-body states are 
unphysical. Formally, the unphysical states do not obey 
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the condition on total fermion parity imposed by the pro¬ 
jector. 

Although not the main focus of this work, it is interest¬ 
ing to repeat the analysis with open boundary conditions. 
More generally, we may consider a lattice with formally 
periodic boundary conditions, but allow for an arbitrary 
number of “missing” bonds with zero bond strength Jif, 
this includes the cases of both open and cylindric bound¬ 
ary conditions. 

A missing a bond, connecting sites i and j, induces two 
dangling gauge Majorana fermions, hf and 6“. These can 
be combined into a canonical fermion, Eq. (161, which 
is decoupled (for zero external field), hence represents a 
zero-energy mode. Occupying this zero mode obviously 
changes the total fermion parity without changing ob¬ 
servable properties of the many-body state. As a result, 
a given Majorana many-body state can always be turned 
from physical to unphysical or vice versa by changing 
the zero-mode occupation. Phrased differently, all matter 
Majorana states in any flux sector are physical if there is 
at least one missing bond which can “absorb” the fermion- 
parity condition. In Appendixj^we demonstrate this for 
a small 2x2 system. A consequence is that the number 
of fermion zero modes of a Kitaev model with missing 
bonds is smaller by one compared to the number of zero 
modes suggested by its Majorana representation. 


IV. SPIN CORRELATIONS AND MAGNETIC 
SUSCEPTBILITY 


the Majorana hopping on the a-bond, representing the 
insertion of the flux pair. A suitable Lehmann represen¬ 
tation of Eq. (211 is in terms of the matter Majorana 
eigenstates of the Hamiltonian T-Lq -I- Ea, denoted by |A): 


= - *5] (MP| c. |A) (A| c, IMP) 

A (22) 

X (5[a; - {Ex - E^)]S(^ij)^Sai3. 

Here, Eq and Ex are the energies of the initial and inter¬ 
mediate states. In the following, the complete sum over 
excited states |A) will be approximately evaluated using 
states with a fixed (small) number of matter excitations 
of Ho + Va', this is a suitable strategy provided that no 
orthogonality catastrophe occurs.l^ 

In order to evaluate the matrix elements {M^lci |A), 
involving eigenstates of both Ho + Va and Ho, we need 
a conversion for the excitation operators. In the follow¬ 
ing we denote the operators for matter eigenmodes in the 
zero-flux and two-flux sectors with d and b, respectively. 
As in Eq. (|^, these are constructed from the matter Ma¬ 
jorana operators according to 

(ai,...,aAr) = ^ [{c^)U + i{cl)V] , 

; (23) 

{h,...,bN) = -[{ci)u' + i{ci)r]. 

Using a Bogoliubov transformation, one can express the 
one kind of operators in terms of the other 


Dynamical spin correlations in the Kitaev model have 
been calculated in Ref. |5] In this section we summarize 
and extend the required formalism. 

Consider the zero-temperature spin correlation func¬ 
tion 

S^^{t) = {0\ar{t)a^{0)\0) (20) 

where |0) is the many-body ground state. Given that 
the fluxes are constants of motion, the correlator can be 
calculated by decomposing the ground state |0) as a di¬ 
rect product of the ground states in the gauge and mat¬ 
ter sector. Specifically, the application of a df operator 
changes the two flux variables which involve the a bond 
emanating from site i. This leads to the dynamical rear¬ 
rangement of matter fermions in the modified gauge field. 
The spin correlator can therefore be expressed purely in 
terms of matter fermions in the ground-st ate flux sector, 
subject to a perturbation Va = —2iJ°‘ciCj^^ 


= (24) 

m 

where X, Y are the transformation matrices 

A* = ^{U'^^U + V'W), 

r (25) 

y* = -{u"^u -V’^V) 

which obey the conditiond^ 

AVI’ -k YY^ = 1, Ar"^ -k YX^ = 0, 

t T . T . i (26) 

A^A-ky^F* = l, X^Y*+Y^X = 0. 

This allows one to rewrite the fermion-free state of the 
two-flux sector, |Ao), in terms of d fermions and the 
fermion-free state in the zero-flux sector, |Mo): 

|Ao) = [AtA]|Mo), (27) 


S^^it) = -t (MP| \M^) 

where Ho is the Majorana hopping Hamiltonian in the 
zero-flux sector and |Mq ) its physical ground state. Site- 
off-diagonal contributions vanish beyond nearest neigh¬ 
bor pairs indicated by {ij)a- Site-diagonal terms are cal¬ 
culated similarly. Hq -k Va and Ho differ in the sign of 


with the overlap |(Mo|Ao)| = i/ldet A|.l^ 

However, as we have pointed out in Section [IH[ in the 
gapless phase the physical states in the flux-free sec¬ 
tor must have an odd number of a fermions. Hence, 
|Mq) = a| |Mo) and Eq = AQ°^-k2e® where and 
are the energies of the ground state and the lowest excita¬ 
tion of Hu in the flux-free sector. Using Eq. (191 we find 
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that I A) must contain an even number of matter fermion 
excitations, see Appendix]^ These multi-particle eigen¬ 
states of "Ho + Va are given by |A) = I "*^ 0)1 

n even. The simplest contribution to (Mg |ci |A) is the 
zero-particle contribution: 


(Mol aiCA,i |Ao) = \/|det A| 


U^o-iUX-^Y)^^ 


(28) 


written for i on the A sublattice. The two-particle con¬ 
tributions can be obtained by straightforward algebra as 


{Mo\aiCA,ib{^b{^ |Ao) = 

V\detX\[u.o Xx.o- 

W. - ^A.o] - 

[X^l-X,,o] - (UX-^Y)^^ (yxn,,,, 

y^^o(Uyn,,^-rMo(UY^hJ- (29) 


Matrix elements for cbj are calculated similarly.!^ 

In contrast, upon ignoring the fermion parity condi¬ 
tion one may start with the fermion-free state |Mo) in 
the zero-flux sector. Then, the spin correlation function 
starts with the one-particle contribution: 


(Mol caJI |Ao) = ^^d^ (30) 


and the energy Eq appearing in Eq. ( [^ is given by Eq^'^ . 

Below we will show results for the dynamic structure 
factor at momentum q = 0 , 

5““(q = 0,a;)=^5““(a.) (31) 

ij 

and the static susceptibility Xiji obtained via the 
Kramers-Kronig relation 

X“/ {u: = 0) = -V [ Au:' , ( 32 ) 

where V denotes the Cauchy principal value. 


V. NUMERICAL RESULTS: CLEAN SYSTEM 

Applying the methodology outlined so far, we now ex¬ 
hibit a few numerical results for the clean Kitaev model, 
obtained via singular-value decomposition of the matrix 
M in Eq. (§. We have treated finite-size systems with 
A 1.2 < 150. Unless noted otherwise, the magnetic cou¬ 
plings are chosen to be isotropic, = J. 


A. Finite-size behavior of the flux gap 

In Fig. we show the finite-size scaling of the energy 
necessary to create a flux pair. Since the physical flux- 
free ground state contains one matter fermion excitation. 



FIG. 2: Flux gap AE of the isotropic Kitaev model as func¬ 
tion of inverse system size, with Li = L 2 = L, periodic bound¬ 
ary conditions, and M = 0. The solid line shows the physical 
result, Eq. ( |33[ ), taking into account the presence of an ex¬ 
cited matter fermion in the flux-free sector. In contrast, the 
dashed line shows the result (341 where both the states in the 
flux-free and two-flux sectors are unphysical. AEp = AE„ is 
realized for L mod 3 = 0 where the Dirac point is an allowed 
wavevector. The arrow indicates the infinite-system resuHP 
AE ^ 0.26J 


whereas the lowest two-flux state does not, the physical 
energy gap is given by 

AEp = E^^'> - El = E^^'> - (4°^ + 2e®) (33) 

where Eq^'^ and E^^'^ are the ground-state energies of Hu 
in the zero-flux and two-flux sectors, respectively, and 
refers to the lowest singular value of M in the flux-free 
sector. Alternatively, one may consider an unphysical 

gap, 

AE, = (4^) + 2ep V 4°^ (34) 

which involves states with incorrect fermion parity in 
both flux sectors. 

As the L = 00 matter fermion spectrum is gapless, we 
have = 0 and thus AEp = AE^, whenever the 

Dirac point is included in the discrete set of momenta. 
For M = 0 this applies to L mod 3 = 0- these data 
points display weak L dependence in Fig.[^ In contrast, 
the data points for L mod 3 yf 0 are influenced by the 
strong L dependence of or \ We note that the 
result in Fig. is qualitatively similar to that in Fig. 4 
of Ref. |T7] where different boundary conditions were em¬ 
ployed. 

Fig. ^ demonstrates that observables calculated for 
physical and unphysical states have rather different 
finite-size behavior; in particular the finite-size conver¬ 
gence appears significantly slower in the physical case. 
Knowing that both AEp and AEu have to converge to 
the same value as L —>■ 00 , one may choose the most 
suitable set of states and boundary conditions for fast 
convergence. 
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B. Ground-state flux sector 

A remark is in order concerning negative values of the 
flux gap for small L, Fig. which imply that the ground 
state is not flux-free. It has been arguecP that a theo¬ 
rem of Liebj^ being concerned with free-particle hopping 
Hamiltonians, guarantees that the ground state of the 
Kitaev model is always in the flux-free sector. This as¬ 
sertion is apparently incorrect, Fig.[^ and the reasons are 
twofold: (i) The theorem of Lieb applies to ground states 
of hopping Hamiltonians, but as established in Ref. |T7| 
and here, the physical ground state of the Kitaev model 
may contain an excited matter fermion which changes the 
energetics (and in particular lowers the energy of the low¬ 
est many-body state in the two-flux sector relative to that 
in the flux-free sector), (ii) Only systems with L 2 = M 
obey the particular periodicity requirement needed for 
Lieb’s theorem to apply. Taken together, the theorem of 
Lieb ensures that the ground state of the Kitaev model 
is in the flux-free sector in the limit of large system size 
(where the restrictions (i) and (ii) become irrelevant), but 
is not decisive for small systems. 


C. Finite-size behavior of the dynamic 
susceptibility 

Fig. 1^ shows the dynamical structure factor calculated 
for two system sizes. Reasonable finite-size convergence 
is apparent,!^ and the results for L = 140 are very close 
to the infinite-system result from Ref. - the latter is 
known to have a gap of size AE/J « 0.26, the flux gap. 

Let us briefly discuss the difference between the phys¬ 
ical and unphysical results. As explained in Section EYl 
the physical flux-free ground state comes with one matter 
fermion excitation, such that (at the isotropic point) the 
excited intermediate states in the two-flux sector have an 
even number of matter fermions. In particular, there is a 
contribution from the zero-fermion intermediate state - 
this produces an isolated S peak in S{uj) at low energies 
(clearly visible in the L = 40 data at w/J « 0.08). The 
rest of the signal comes from two-fermion intermediate 
states; higher excited states are ignored in our calcula¬ 
tion because they only carry spectral weight of about 
2.5%.^ In contrast, the unphysical signal is obtained by 
starting from a fermion-free ground state in the flux-free 
sector. Then, the signal at the isotropic point arises from 
single-fermion intermediate states, and the low-energy S 
peak is absent. 

Remarkably, the differences between the physical and 
unphysical signal diminish with increasing system size, in 
accordance with the general argument from Section [HIB| 
Here, the reason for this can be understood in detail: 
Although the two-fermion intermediate states |A) in the 
physical case can have two arbitrary fermions excited, 
the matrix element (A| Ci \Mq) will only be sizeable if 
one of the fermions is the lowest-energy one, simply 
to match the lowest-energy fermion occupied in \Mq). 
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FIG. 3: Dynamic structure factor for the isotropic Kitaev 
model, calculated from Eq. (311 for systems with Li ,2 = 40 
and a broadening of S/J = 0.04 (top) and Li ,2 = 140 and 
5/J = 0.02 (bottom), both with M = 0. The “physical” 
(solid) result takes into account the presence of a matter 
fermion in the ground state; it consists of two-particle con¬ 
tributions, Eq. (291, and an isolated low-energy peak corre¬ 
sponding to the zero-particle contribution, Eq. (28l. In con¬ 
trast, the “unphysical” result (dashed) contains one-particle 
contributions, Eq. (301, only. The exact resullPfor L = 00 is 
shown for comparison. 


All other matrix elements are suppressed at least with 
which effectively reduces the two-particle contin¬ 
uum to the single-particle continuum of the unphysical 
case. Similarly, the matrix element (Ao|ci|MQ), deter¬ 
mining the weight on the low-energy 6 peak in the phys¬ 
ical response, scales as Hence, the dynamical 

structure factor in the thermodynamic limit is indepen¬ 
dent of the ground-state parity tt.I^ 


VI. BOND DISORDER: GENERAL 
CONSIDERATIONS 


Before showing numerical results for the Kitaev model 
with bond disorder, we quickly summarize a few general 
aspects, some of which have been discussed in Refs. 16171 

133 



























Provided that the ground state in the presence of dis¬ 
order remains in the flux-free sector, the low-energy be¬ 
havior in the presence of bond disorder is equivalent 
to that of Dirac fermions with random hopping on the 
honeycomb lattice. This is a special case of a bipar¬ 
tite random-hopping problem, belonging to the symme¬ 
try class BDI in the Altland-Zirnbauer classiflcation.l^ 
The single-particle properties of such systems have been 
analyzed using various techniquesP^HMl ^11 single-particle 
states at non-zero energies are exponentially localized, 
and the result ing density of states at low energies follows 
the forn J'^^P^ J 


p{uj) cx — exp (—c\ Inwl^^^') 

UJ \ / 


(35) 


with X = 3/2. This immediately implies a corresponding 
singular behavior for the specific-heat coefficient CjT. 
However, the asymptotic form (35) is only realized be¬ 


low an extremely small energy scale which depends on 
the disorder strengtlP^ and is typically not accessible in 
numerical simulations. 

In the application to the Kitaev model, two further as¬ 
pects are important: (i) For strong disorder, the ground 
state may not be located in the flux-free sector - this 
will be discussed in Section VIII (ii) Even if the ground 
state is in the flux-free sector, the flux gap A may be¬ 
come small, and many-body states in excited flux sectors 
become important for temperatures T ^ A. 


VII. NUMERICAL RESULTS: DISORDERED 
SYSTEM 



- 0.1 0 0.1 0.2 0.3 0.4 0.5 


AE[J| 



FIG. 4: Distribution of the local flux gap for the isotropic 
Kitaev model with bond disorder, calculated for L = 40 (top) 
and L = 80 (bottom) and two different values of the disorder 
strength A“. Shown are the results for both the “physical” 
(closed symbols) and the “unphysical” (open symbols) gap, 
calculated according to Eqs. (|33[) and (|34[), respectively. 


A. Flux gap 


B. Static susceptibility 


Figure]^ shows histograms of the local flux gap, AEij, 
for Kitaev models with bond disorder. This local gap is 
defined as in Eqs. ([33|) and (341, with the specific two-flux 


state obtained by flipping the{ij) bond. We note that the 
selection rules for physical states continue to apply for the 
moderate disorder considered here. Comparing the L — 
40 and L = 80 data, strong finite-size effects are apparent 
which are inherited from the disorder-free situation, see 
the results in Fig.[^ The following discussion thus mainly 
applies to the L = 80 data. 

For weak disorder, A“/J = 0.1, the gap distribu¬ 
tion is essentially symmetric, with a relative width which 
roughly matches that of the coupling-constant distribu¬ 
tion. For strong disorder, A“/J = 0.5, the gap dis¬ 
tribution widens and becomes slightly asymmetric. Its 
mean value is shifted downwards relative to the clean case 
(there AEp/J = 0.173 and AE^/J = 0.345 for L = 80). 
Furthermore, cases with AE < 0 appear, i.e., the ground 
state is not in the flux-free sector. The significance of 
this finding will be discussed in Section |VIII[ 


With an eye towards nuclear-magnetic-resonance ex¬ 
periments, we consider the local susceptibility 

XNMR(j) = (36) 

3 

which is proportional to the resonance frequency in NMR 
experiments. We recall that, in the Kitaev model, Xij = 0 
beyond nearest-neighbor distance, i.e, there are only on¬ 
site and nearest-neighbor contributions to xnmr- 

Results for the distribution of Xnmr(*) are displayed 
in Fig. 1^ While weak disorder again produces an essen¬ 
tially symmetric distribution with a relative width cor¬ 
responding to that of the coupling-constant distribution, 
strong disorder produces a distinctly asymmetric shape 
with a tail at large values of y. The reason is in the 
strong fluctuations of the flux gap. Fig. considering 
that X (X 1 /AE. We note that in evaluating x we have as¬ 
sumed the ground state to be flux-free, and consequently 
have discarded the rare events with AE < 0. 
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FIG. 6: Dynamic structure factor as in Fig. but now for 
the Kitaev model of size L\ = L 2 = 80 with box-type bond 
disorder of strength A“/J = 0.5. The artificial broadening is 
smaller than in Fig. d/J = 0.01. The clean-system resullP 
is shown for comparison. 


tween the physical and unphysical results at L = 80 in 

Fig.i 


VIII. TRANSITION OUT OF FLUX-FREE 
STATE 


FIG. 5: Distribution of the local (NMR) susceptibility, 

Eq. (361, for the isotropic Kitaev model with bond disorder, 
calculated for L = 40 (top) and L = 80 (bottom) and two 
values of the disorder strength A“. As before, “physical” (“un¬ 
physical”) represent the results obtained for one (zero) matter 
fermions in the flux-free ground state. 


Interestingly, and in striking contrast to the results for 
the flux gap in Fig. we And that the physical and un¬ 
physical results for the y distribution are almost identical 
at L = 80. The explanation is similar to that given in 
Section |V C| Although the physical and unphysical cases 
have contributions to y with rather different excitation 
energies, the corresponding matrix elements are small for 
large L. For instance, the zero-particle contribution to 
the physical susceptibility, with the excitation energy be¬ 
ing the flux gap AAp according to Eq. (331, has a weight 
scaling as N~^. 


C. Dynamic susceptibility 

As a further example, we plot the dynamic structure 
factor in the presence of bond disorder in Fig.|^ As disor¬ 
der tends to smear the flux gap, the gap in the structure 
factor is filled. This is accompanied by a shift of weight to 
lower energies, as expected from Fig.|^ Disorder-induced 
changes at higher energies are minimal. Consistent with 
the above discussion, there is essentially no difference be- 


Our numerical results show that, with increasing bond 
disorder, the ground state of a finite-size Kitaev model 
is no longer located in the flux-free sector. Instead, the 
ground state displays a finite flux density, where fluxes 
occur in the system at random positions which depend 
on the disorder realization. We note that such a state is 
trivially realized for box disorder with A“/J > 1, as this 
implies the existence of bonds with flipped sign which can 
be compensated by placing flux pairs adjacent to these 
bonds (equivalent to choosing u = — 1 on the respective 
bonds). More interesting is the possible occurrence of 
such a random-flux state for A“/J < 1 where all bond 
strengths are positive. Notably, the numerics also indi¬ 
cates that the tendency towards ground-state fluxes di¬ 
minishes with increasing L (see e.g. Fig. |^, such that 
definite conclusions about the thermodynamic limit can¬ 
not be drawn. 

However, we are able to provide a general argument in 
favor of a non-trivial transition to a random-flux state 
which applies to the thermodynamic limit. A key ingre¬ 
dient is the observation of Refs. iftti that a vacancy site 
gains a finite amount of energy by binding a flux. Con¬ 
sider now the more general situation where a single defect 
site is surrounded by three bonds of strength J' and em¬ 
bedded in an otherwise homogeneous Kitaev model with 
couplings J. While J' = 0 corresponds to the vacancy 
case, this site will also bind a flux for finite small J'. 
This is shown in Fig. For all 0 < J' < Jmin with 
Jmi-a/J ~ 0.04 the energy of the state with a flux bound 
in one of the three plaquettes adjacent to the defect is 
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FIG. 7: Flux energy Efiux{J') for an isotropic Kitaev model 
with a single defect site which has three weak bonds of 
strength J' to its neighbors, influx < 0 implies that the de¬ 
fect binds a flux; influx(0)/J = —0.027 is the known result for 
a vacancy from Ref. |6l i5flux( ^0 been calculated as the 
energy difference between a two-flux state, with one flux in 
one of the three defect plaquettes and one flrrx at maximum 
distance away from the defect, and the flux-free state, and 
the energy of an isolated flux (for the same L) has been sub¬ 
tracted. The inset shows the Hnite-size scaling of influx, the 
error bars in the main panel arise from uncertainties in the 
1/ —>■ oo extrapolation. The dashed line is a linear fit. 


lower than that of the flux-free state. More generally, a 
defect site surrounded by three bonds of a strength in the 
interval [0, J„iin] will bind a flux. 

Now, for box disorder with strength A“ the mini¬ 
mum coupling strength is J — A“, thus that for any 
J > J — A“ there is a finite probability to And local con¬ 
figurations which have (i) three bonds emanating from 
one site with strength smaller than J and (ii) all sur¬ 
rounding bond strengths arbitrarily close to J. This 
is exactly the condition for locally binding a flux, pro¬ 
vided that J < Jmin- We conclude that a random- 
flux state must be realized for disorder strengths with 
A“ > J — Jmin- This proves the existence of a transi¬ 
tion - from zero flux to random flux - somewhere in the 
interval 0 < A“/J < 1 — Jmin/J ~ 0.96. 


IX. SUMMARY 

Our study of Kitaev’s honeycomb model with bond 
disorder has lead to twofold results: On the one hand, 
we have dealt with the selection of physical states in 
the Majorana representation. Extending earlier work, 
we have shown that the ground state of the gapless Ki¬ 
taev model with periodic boundary conditions generically 
contains one matter fermion excitation. This causes sig¬ 
nificant finite-size effects for observables, as illustrated 
for the flux gap. We have also discussed the difference in 
state selection between the cases with periodic and open 
boundary conditions. Obviously, this state selection is of 


relevance for all numerical studies of Kitaev models using 
Majorana fermions. It will be interesting to extend this 
analysis to other tricoordinated lattices where the Kitaev 
model can also be solved exactly; work in this direction 
is in progress. 

On the other hand, we have numerically determined 
the static and dynamic spin susceptibility in the pres¬ 
ence of bond disorder. In particular, we have calculated 
the distribution of local susceptibilities which determines 
the NMR lineshape. For large disorder, we predicted a 
transition to a random-flux state. A detailed study of 
this transition is left for future work. 
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Appendix A: Parity of matter fermion excitations 
1. Gapless phase 

The purpose of this appendix is to prove that all flux- 
free physical states in the gapless phase of a translation- 
invariant Kitaev model with periodic boundary condi¬ 
tions contain an odd number of dm fermion excitations. 
This supersedes the results of Ref. |T71 but is consistent 
with their Fig. 3. 

The proof is based on insights from Ref. [IT] which we 
lay out first. The flux-free sector is characterized by all 
Uij = 1. Then the eigenmodes of "Hu are diagonal in 
momentum spaceP 

■Hu = XI - 1) (Al) 

q 

with /(q) = -|- -|- J«. The spectrum 

|/(q)| is gapped if or permutations, and 

gapless otherwise. The reciprocal lattice is defined by 
the vectors hi^ 2 , see Fig. For any finite lattice the 
Brillouin zone is reduced to a finite set of wavevectors 
q, which can be partitioned into three sets H and H±. 
We assign q G H if ±q are equivalent (up to reciprocal 
lattice vectors); there are at most four wavevectors in fl, 
namely 0, bi/2,b2/2, and (bi -|-b2)/2. The remaining 
q are partitioned such that ±q belong to two distinct 
sets . One can then derive the explicit formula for the 
determinant of the transformation matrixETl 

det(g“) =-l'>'+^', (A2) 
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L\ L 2 ^ 

7 

{LiL2f LiM 

(-ir 

+ + + 

+ - - 

+ + - 

+ - + 

- + + 

- + - 

- - -f 

1 

1 

0 

0 

0 

0 

0 

0 

+ + 

+ + 

+ + 

+ + 

+ + 

+ 

-f 

-1 

-1 

-1 

-1 

-1 

-1 

-1 

-1 


TABLE I: This table shows (-l)®det(g“) = (-1)'' for the 
gapless phase in relation to the boundary conditions I/i, 2 , M 
(where + and — refer to even and odd values, respectively) 
and the resulting 7 , see text. 


valid for the flux-free sector, where N = L 1 L 2 , and 7 
is the number of reciprocal vectors q S with /(q) < 
0. Together with the geometric factor (14), we can now 
rewrite 

(-l)®det(Q") = (_i) 7 +ii+i 2 +i?ii+iiM-M 2 ^ (- 1 )'". 

(A3) 

Although 7 depends in a non-trivial way on the bound¬ 
ary conditions Li ^2 and M as well as on the couplings 
Jx,y,z, we can calculate it for any given choice of Li, 2 , 
M. Since hiBj = 27r5y, it is easy to see that in the gap¬ 
less phase only / — jy ig less then 0 . 

Therefore 7 = 1 if (bi -|- b2)/2 G il and 7 = 0 otherwise. 
The allowed q vectors are determined by the conditions 


g^qilbi ^ 
giq(L2b2-|-Mbi) _ ^ 


(A4) 

(A5) 


with q = qi^i -|- g 2 b 2 - Therefore 7 = 1 if Li = 2ni 
and L 2 + M = 2n2 (ni .2 G Z). Enumerating all eight 
combinations of parities of Li ^2 and M yields the results 


in table p 
equation 


showing that —1^ = —1 in all cases. Using 
1 ^ this implies that, in the flux-free case where 


= 0, the physical Majorana states must have an odd 

number of matter fermion excitations, tt = (— 1 )'^“ = — 1 . 

From this result one can further deduce that Na for 
states in the two-flux sector, at and near the isotropic 
point, is even. This flux sector has = 1, and Eq. (19) 

(— 1 )® det(( 5 “)(— 1 )^^ (— 1 )^“ = 1 implies that Na must 
be even as long as the signs of det(Q“) in the zero-flux 
and two-flux sectors are identical. The latter applies near 
the isotropic point, but not in the entire gapless phase . 1 ^ 



J„,Jy [JJ 


FIG. 8 : Lower half of the many-body spectrum of an 

anisotropic 2x2 Kitaev model with as fnnc- 

tion of , with the system geometry shown in the inset. 

Lines: Eigenenergies obtained by exact diagonalization of the 
spin Hamiltonian. Symbols: Eigenenergies of the Majorana 
Hamiltonian in the flux-free sector, Uij = 1. Na is the num¬ 
ber of matter fermion excitations. At (and near) the isotropic 
point, Na = 0,2 states are unphysical (red) while Na = 1 
states are physical (blue). The vertical dashed line indicates 
the boundary between the gapped and gapless phases.® 


'y{Li, L 2 , M) is different from that in the gapless phase, 
and 7 can now take all values from 0 to 3. As a result, 
a unique conclusion similar to the gapless phase cannot 
be reached. Moreover, the small flux gap in combina¬ 
tion with the large fermionic gap can lead to the physical 
ground state having excited flux pairs but no fermions, 
see also Fig. 5 of Ref. [T71 


Appendix B: Spectrum for Li = L 2 = 2 

In this appendix we verify the analysis in Section |III| 
by comparing the eigenenergies of "Hki obtained by exact 
diagonalization of the spin Hamiltonian, with the ener¬ 
gies of the many-body Majorana states, both physical 
and unphysical. 

We choose a small system with Li = L 2 = 2 and M = 
0. Here the Dirac point does not belong to the discrete 
partitioning of the Brillouin zone, such that all excitation 
energies of matter fermions, Cm, are non-zero. 


2. Gapped phase 

Although a similar analysis may be performed for the 
gapped phase of the Kitaev model, it turns out that the 
different parity combinations of L 12 and M come with 
different signs for (—1)^. In particular, the dependence 


1. Periodic boundary conditions and varying 
anisotropy 

To illustrate the unphysical character of the zero-flux 
fermion-free state, we show in Fig.j^the many-body Ma- 
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FIG. 9: Same as Fig.[^ but for the four-flux sector with Wi = 
IF 2 = — 1. The bonds with Uij = — 1 are shown in light (red) 
color in the inset. Here, Na = 0, 2 states are physical (blue) 
while Na = 1 states are unphysical (red) near the isotropic 
point. 

jorana energies in the zero-flux sector, together with all 
2® = 256 eigenenergies of the spin Hamiltonian, for vary¬ 
ing spin anisotropy. 

In the entire gapless phase, 1/2 < < 1, the 



Jo [J\ 


FIG. 10: Same as Fig. but now for an isotropic model 
where a single bond has a different exchange strength Jo 7 ^ J. 
Full (open) symbols correspond to the Majorana eigenener¬ 
gies in the sectors with zero flux (two fluxes, with a flux pair 
adjacent to the Jq bond), respectively. As before, blue (red) 
symbols denote physical (unphysical) states. 


Majorana states with even number Na of matter fermion 
excitations do not correspond to any of the physical 
states, whereas the Majorana states with odd Na match 
the physical spectrum. Interestingly, this behavior is re¬ 
versed in the gapped phase, 0 < < 1 / 2 , where 

now the states with even Na are physical. 


We have repeated this analysis in all flux sectors. As an 
example, we show the flux sector containing the ground 
state, here with fluxes through all plaquettes, in Fig. 
The physical states in this sector have an even number 
of excited matter fermions in both phases. 


Interestingly, in the three flux sectors without plaque- 
tte fluxes but with a flux through at least one of the 
torus holes, i.e., Wi = — 1 , W 2 = 1, Wi = 1 , W 2 = — 1 , 
and Wi = W 2 = — 1, the even-iVa states are found to be 
physical. 


2. Varying a single bond 


To underline the arguments concerning missing bonds 
and open boundary conditions in Section |III C| we now 
consider an isotropic Li = L 2 = 2 system where we vary 
the exchange strength Jq on one bond keeping the other 
couplings fixed at J. Fig.[T0]shows the Majorana energies 
both in the zero-flux and two-flux sectors, in the latter 
case with the flux pair located adjacent to the Jq bond, 
together with the exact spectrum. 


For any non-zero Jq, the states with odd (even) Na are 
physical in the zero-flux (two-flux) sector, respectively, 
consistent with our reasoning above. However, for Jq = 
0, all matter Majorana states become physical: This is 
a consequence of the zero mode constructed from gauge 
Majorana fermions in the presence of a missing bond, see 
Section [III C[ Consistent with this, the energy difference 
between the zero-flux and two-flux states vanishes as the 
flux pair has no observable impact if it surrounds the 
Jo = 0 bond. 
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